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INTRODUCTION 


This  report  is  a brief  description  of  the  main  aspects  of  a cqmputer  simulation 
model  designed  to  mimic  air/ground  communications  in  an  air  traffic  control 
(ATC)  sector.  The  model  was  developed  through  a joint  effort  of  the  National 
Aviation  Facilities  Experimental  Center  (NAFEC)  and  Princeton  University  to 
apply  fast  time  simulation  techniques  and  methods  of  time  series  analysis  to 
ATC  problems  (references  1 to  4).  The  model  has  been  validated  with  field 
data  obtained  in  1969  from  the  New  York  Control  Area  and  additional  informa- 
tion collected  from  the  Houston  Control  Area  in  1971.  The  model  is  presently 
available  in  a software  package  adapted  for  use  in  computer  facilities  at 
NAFEC.  It  can  be  used  to  simulate  the  following  nine  sector  functions: 

HI  High-altitude  enroute  (radar  controllers) 

LE  = Low-altitude  enroute  (radar  controllers) 

LT  ■ Low-altitude  transitional  (radar  controllers) 

GN  * Ground  control  (tower  controllers) 

LC  ■ Local  control  (tower  controllers) 

LG  - Local  ground  control  (tower  controllers) 

DP  - Radar  departure  control  (IFR  (instrument  flight  rules)  Room  radar 
controllers) 

AD  ■ Radar  arrival/departure  control  (IFR  Room  radar  controllers) 

AR  ■ Radar  arrival  control  (IFR  Room  radar  controllers) 

In  the  following  paragraphs,  we  describe  the  structure  of  the  model,  input 
variables,  output  variables,  and  some  possible  applications.  More  detailed 
descriptions  can  be  found  in  references  2 and  3. 


DISCUSSION 


STRUCTURE  OF  COMMUNICATIONS. 

Before  proceeding  to  a discussion  of  the  model,  some  formulation  of  the 
structure  of  air/ground  communications  is  in  order.  An  air/ground  conversa- 
tion normally  consists  of  several  transmissions  (TR's)  alternately  initiated 
by  pilot  and  controller.  In  keeping  with  earlier  work,  a whole  conversation 
is  referred  to  as  a communication  transaction  (CT),  so  that  a CT  consists  of 
one  or  more  TR's.  The  reader  is  advised  to  eep  the  distinction  between  TR's 
and  CT's  in  mind. 

Obviously,  while  an  aircraft  is  in  sector,  there  are  several  CT's  between  con- 
troller and  pilot  as  shown  In  figure  1.  From  the  standpoint  of  communications, 
we  adopt  the  attitude  that  an  aircraft  arrives  at  a sector  with  the  beginning 
of  the  first  CT  with  the  ground  and  leaves  the  sector  at  the  end  of  the  final 
CT.  Between  these  points  in  time,  the  communication  between  a single  aircraft 
and  ground  usually  consists  of  several  CT's  and  gaps  between  CT's.  We  refer  to 
these  gaps  as  intercommunication  gaps.  Referring  to  figure  1,  it  is  apparent 
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that  the  number  of  intercommunication  gaps  is  one  less  than  the  number  of  CT's, 
so  that  specification  of  either  automatically  determines  the  other. 

Figure  1 is  essentially  a picture  of  air/ground  communications  from  the  point 
of  view  of  a pilot.  Since  there  are  usually  many  aircraft  in  sector  simulta- 
neously, the  communications  picture  from  the  point  of  view  of  the  controller 
is  just  a superposition  of  several  diagrams  like  that  of  figure  1.  This  is 
demonstrated  in  figure  2 for  the  case  where  the  maximum  number  of  aircraft  in 
sector  between  time  1 and  2 is  two.  Obviously,  the  gaps  between  CT's  from  the 
viewpoint  of  the  controller  are  normally  shorter  than  those  experienced  by  the 
pilot  of  a single  aircraft.  In  order  to  distinguish  the  former  from  the  latter, 
we  refer  to  the  gaps  experienced  by  the  controller  as  intertransaction  gaps. 
Thus,  the  pilot  experiences  intercommunication  gaps,  and  the  controller 
experiences  intertransaction  gaps. 

STRUCTURE  OF  THF.  PROGRAM. 

The  simulation  model  exists  as  a computer  program  written  in  the  GPSS  V and 
FORTRAN  IV  languages.  It  is  designed  to  mimic  second-by-second  behavior  of 
sector  air/ground  communications  over  periods  of  time  in  the  order  of  hours. 

A flow  chart  of  the  model  is  displayed  in  figure  3.  As  indicated  by  the  blocks 
of  the  chart,  the  model  performs  nine  basic  operations.  These  are  described 
below. 

Aircraft  arrivals  (block  1)  at  the  sector  under  investigation  correspond  to  a 
sample  function  from  a Poisson  process.  Stated  another  way,  interarrival 
times  are  modeled  as  independent  exponential  variates  with  a common  average 
(AMEAN)  expressed  in  seconds.  AMEAN  is  an  input  variable  that  specifies  the 
average  rate,  in  this  case,  1/AMEAN  aircraft  per  second,  with  which  aircraft 
enter  the  sector. 

When  an  aircraft  arrives  at  the  sector,  it  is  assigned  a number  of  CT's 
(block  2)  by  means  of  a random  sample  drawn  from  a negative  binomial  distribu- 
tion with  shifted  origin  and  parameters  P and  K,  i.e., 

K+r-2  pK(l-P)r-1;  r-1,2,3, 

K-l 

For  example,  P and  K might  be  estimated  by  the  method  of  moments  from  histori- 
cal data  collected  from  a sector  or  group  of  sectors  of  the  type  being  studied. 
In  addition  to  the  number  of  CT's,  an  incoming  aircraft  is  assigned  a mean 
intercommunication  gap  length  (MGAP)  (block  3).  The  natural  logarithm  of 
MGAP  is  generated  from  a normal  distribution.  The  mean  (XM)  of  this  distri- 
bution is  given  by, 

XM  - A1  + A2XN 

where  N is  the  number  of  gaps,  one  less  than  the  number  of  CT's,  and  A1  and  A2 
are  regression  coefficients  determined  from  historical  records  of  mean  gap 
length  and  number  of  gaps  obtained  for  many  aircraft  passing  through  sectors 
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CT  (DURATION  VARIES) 
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FINAL  CT  ON  LEAVING  SECTOR 
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FIGURE  1.  COMMUNICATIONS  FROM  PILOT  POINT  OF  VIEW 


SINGLE  AIRCRAFT  COMMUNICATIONS  PERFORMANCE 
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COMPOSITE  COMMUNICATIONS  PERFORMANCE 


TIME  IN  SECONDS 


FIGURE  2.  COMMUNICATIONS  FROM  CONTROLLER  POINT  OF  VIEW 
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FIGURE  3.  BLOCK  DIAGRAM  OF  SIMULATION  MODEL 


4 


of  the  type  being  studied.  The  standard  deviation  (SD)  is  determined  from  the 
relationship. 


CU-XM 

2.5758  if  A 2 >0 


SD  = 


I XM-CL 
' 2.5758 


if  A2  < 0 


where  CU  and  CL  are  upper  and  lower  bounds,  respectively,  of  the  logarithm  of 
recorded  mean  gap  lengths.  These  formulas  are  the  result  of  examinations  of 
historical  records  that  indicate  the  existence  of  patterns  in  the  data  of  the 
type  shown  in  figures  4 and  5.  The  constant  2.5758  was  chosen  to  insure  that 
99  percent  of  the  probability  mass  of  the  normal  distribution  used  to  generate 
the  logarithm  of  the  mean  gap  length  lies  between  2 XM  - CU  and  CU  in  the  case 
where  A2  >0,  and  between  CL  and  2 XM  - CL  when  A2  <0.  In  either  case,  whenever 
sampling  from  this  distribution  results  in  a negative  number,  the  mean  gap 
length  is  arbitararily  set  equal  to  1 second. 

The  first  CT  assigned  to  an  aircraft  commences  as  the  aircraft  enters  the  sec- 
tor provided  that  the  channel  (block  6)  is  available,  i.e.,  a conversation 
between  another  aircraft  and  ground  is  not  taking  place.  Otherwise,  the  CT 
enters  a queue  (block  5)  on  a first-come,  first-serve  basis.  As  each  CT  is 
completed,  the  model  ascertains  whether  or  not  all  assigned  CT’s  have  taken 
place  (block  7).  If  not,  then  an  intercommunication  gap  length  is  randomly 
selected  from  an  exponential  distribution  with  mean  MGAP  seconds  (block  8). 

If  it  is  determined  that  the  CT  just  completed  is  indeed  the  last  of  the 
assigned  CT’s,  then  the  aircraft  leaves  the  sector  (block  9).  At  any  rate, 
the  end  time  of  the  most  recently  completed  CT  and  the  length  of  the  sub- 
sequent intercommunication  gap  determines  the  start  time  of  the  next  CT.  Of 
course,  in  the  event  that  the  channel  is  busy  when  the  CT  is  scheduled  to 
take  place,  then  the  CT  enters  a queue  and  waits  in  turn  for  transmission  ser- 
vice. When  the  start  time  of  a CT  is  established,  the  corresponding  CT  length 
is  determined  in  a two-step  process  (block  4).  The  first  step  of  the  process 
is  to  establish  the  number  of  TR's  involved  in  the  CT.  This  is  accomplished 
by  a random  sample  drawn  from  an  empirical  distribution  based  upon  data  col- 
lected from  one  or  more  sectors  of  the  type  being  investigated.  Thereafter, 
the  length  of  each  TR  is  obtained  as  a random  sample  drawn  from  a gamma  dis- 
tribution with  parameters  a and  A,  i.e., 
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FIGURE  5. 
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The  parameter  values  are  obtained  from  two  master  equations  (1)  and  (11). 
These  are 


(i)  a 

. X =3.70  - 

56. 88/ AMEAN 

(General) 

=2.97  - 

15 . 12/AMEAN 

(LC) 

(ii) 

( a -1)A  =2.0 

(General) 

( a -1)X  =1.7 

(GN) 

( a -1 . 13)A  =1 .72 

(DP) 

( a -1 . 5)A  =1.0 

(AD) 

where,  as  already  indicated,  AMEAN  Is  the  mean  aircraft  interarrival  time  in 
seconds.  Thus,  if  the  LC  sector  function  is  to  be  simulated,  then  appropriate 
values  for  a and  Xare  obtained  from  the  two  equations, 

a • A =2.97  - 15 .12 /AMEAN 

(a  -1)  X =2.0 

A listing  of  the  program  available  at  NAFEC  is  provided  in  the  appendix. 

This  listing  together  with  the  description  of  the  program  structure  presented 
here  is  merely  intended  to  provide  the  reader  with  some  general  idea  of  the 
existing  simulation  capability.  Further  details  including  possible  modifica- 
tions, variations  in  the  manner  in  which  inputs  can  be  supplied  to  the  model, 
justification  of  model  formulation,  etc.,  can  be  obtained  from  the  cited 
references . 

INPUT  VARIABLES . 

As  described,  the  model  is  characterized  by  10  input  variables,  namely,  AMEAN, 
P,  K,  Al,  A2,  CU,  CL,  a,  X,  and  an  empirical  distribution  for  the  number  of 
TR's  contained  in  one  CT.  AMEAN  prescribes  the  mean  interarrival  time  between 
successive  aircraft  entering  the  sector.  P and  K determine  the  distribution 
of  the  number  of  CT's,  or  equivalently,  the  number  of  intercommunication  gaps, 
experienced  by  an  aircraft  as  it  passes  through  the  sector.  The  four  param- 
eters, Al,  A2,  CU,  and  CL,  determine  the  distribution  of  gap  length.  Finally, 
a,  X,  and  the  empirical  distribution  of  the  number  of  TR's  in  a CT  specify 
the  distribution  of  CT  length. 

In  an  application  of  the  model,  all  input  parameters  except  AMEAN  might  be 
assigned  fixed  values  to  represent  a particular  sector  in  some  center  such 
as  Houston  or  New  York.  Then  several  simulations  could  be  performed  cor- 
responding to  decreasing  values  of  AMEAN  to  ascertain  the  effect  of  increasing 
aircraft  arrival  rate  on  sector  communications.  By  the  same  token,  AMEAN 
could  be  held  constant  and  the  effect  of  some  other  input  parameter  on  sector 
communications  observed. 
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RESPONSE  VARIABLES. 


Three  basic  response  variables  are  observed  during  each  second  of  simulated 
time.  They  are  the  sector  aircraft  loading  (i.e.,  the  number  of  aircraft  that 
have  been  handed  off  to  the  sector,  but  have  not  yet  been  handed  off  bjr  the 
sector),  the  state  of  the  air/ground  channel  (i.e.,  busy  or  idle),  and  the 
queuing  state  of  each  aircraft  in  the  sector  (i.e.,  in  queue  or  otherwise). 

An  aircraft  is  defined  to  be  in  queue  if  either  the  controller  or  pilot  desires 
to  converse  with  the  other,  but  is  prevented  from  doing  so  by  virtue  of  the 
fact  that  a conversation  is  presently  taking  place  between  another  aircraft 
and  ground.  Thus,  insofar  as  the  computer  simulation  is  concerned,  the  third 
response  variable  measures  the  lag  between  the  instant  that  the  first  TR  of  a 
CT  is  scheduled  to  occupy  the  channel  and  the  time  that  it  is  actually  carried 
by  the  channel.  In  this  sense,  the  queuing  state  of  the  channel  represents 
delay  in  the  transfer  of  information  between  air  and  ground. 

The  next  few  paragraphs  discuss  each  basic  response  variable  in  greater  detail 
in  the  context  of  a specific  simulation  run.  Values  selected  for  input  para- 
meters in  this  example  are  as  follows: 

AMEAN=86  seconds  P=0.495  K=3.88 

Al=4 . 336  A2=0 .032  CU=6.0 

CL=3 . 1 

The  example  experiment  was  designed  to  simulate  New  York  LC  sector  510.  As  a 
result,  a and  A were  determined  from  the  second  of  master  equations  (i)  and  the 
first  of  master  equations  (ii).  Moreover,  the  empirical  distribution  of  the 
number  of  TR's  in  a CT  was  determined  from  historical  records  obtained  from 
sector  510.  The  simulation  was  allowed  to  run  for  1 hour  of  simulated  time 
prior  to  the  accummulation  of  any  data  in  order  to  dissipate  the  influence 
of  transient  phenomena  generated  by  boundary  conditions  that  exist  at  the 
beginning  of  the  experiment.  Thereafter,  data  were  gathered  for  2 hours  of 
simulated  time.  Consequently,  results  obtained  for  the  2-hour  observation 
period  can  be  viewed  as  representative  of  sample  functions  drawn  from  station- 
ary random  phenomena.  Of  course,  in  the  event  that  a transient  effect  persists 
or  that  values  assigned  to  input  parameters  result  in  an  explosive  situation, 
then  the  underlying  stochastic  processes  are  far  from  stationary.  However, 
such  circumstances  are  usually  reflected  by  very  definite  trends  in  the  time 
series  generated  by  one  or  more  of  the  basic  response  variables  during  the 
2-hour  observation  period. 

CHANNEL  UTILIZATION  AND  AIRCRAFT  LOADING. 


From  observation  of  the  first  response  variable,  it  is  possible  to  compute  the 
number  of  aircraft  in  sector  per  second  averaged  over  each  minute  of  simulated 
time.  The  average  over  the  t th  minute  of  simulated  time  is  represented  by 
nt,  and  there  are  120  such  averages;  namely,  n^  through  n^O*  These  are 
illustrated  for  LC  sector  510  in  figure  6.  The  number  of  aircraft  in  sector 
averaged  over  2 hours  of  simulated  time  is  lust  the  arithmetic  mean  of 
ni  through  n^O'  For  example,  the  2-hour  average  corresponding  to  figure  6 
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FIGURE  6.  SIXTY-SECOND  AVERAGE  AIRCRAFT  LOADINGS 
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is  5.983  aircraft  per  second,  and  this  is  listed  in  table  1 along  with  the 
maximum  number  of  aircraft  in  sector  during  any  second  of  simulated  time  and 
the  total  number  of  aircraft  handled  by  the  sector  during  the  2-hour  sample 
period.  Along  the  same  lines,  second-by-second  observations  of  the  channel 
state,  busy  or  idle,  can  be  used  to  determine  the  fraction  of  the  2-hour  sample 
period  during  which  the  channel  is  busy  as  well  as  the  fraction  of  the  t th 
minute  of  simulated  time  during  which  the  channel  is  utilized.  The  latter 
fraction  is  denoted  by  cj;  and  is  called  the  channel  utilization;  the  former  is 
called  the  average  channel  utilization.  The  120  values  of  channel  utilization 
are  shown  in  figure  7.  The  average  channel  utilization  is  displayed  in  table  1 
together  with  the  total  number  of  air/ground  CT's  generated  during  the  2-hour 
period,  and  the  average  length  of  these  CT's. 

COMMUNICATIONS  DELAY . 


We  now  turn  to  the  last  of  the  three  basic  response  variables;  namely,  the 
aircraft  queuing  state.  As  an  aircraft  enters  the  sector,  the  simulation 
software  assigns  the  number  of  CT's  that  are  to  take  place  while  the  aircraft 
is  subject  to  sector  control  and  the  start  time  of  the  initial  CT.  Thereafter, 
upon  completion  of  a CT,  the  software  schedules  the  start  time  of  the  next 
CT,  if  any.  In  the  event  that  two  or  more  CT's  vie  for  simultaneous  use  of 
j the  channel,  a queue  of  CT's  is  established.  As  a result,  the  scheduled  start 

time  of  a CT  and  i.he  actual  start  time  need  not  coincide.  As  already  Indicated, 
an  aircraft  is  considered  to  be  in  queue  whenever  a scheduled  CT  between  that 
aircraft  and  ground  is  in  queue.  Thus,  second-by-second  observations  of  the 
aircraft  queuing  state  can  be  used  to  compute  the  delay  between  the  instant 
that  the  CT  would  take  place  in  the  presence  of  an  unlimited  number  of  channels 
and  the  time  that  the  air/ground  channel  does  in  fact  become  available.  Obviously, 
the  delay  is  dependent  upon  the  queue  discipline  incorporated  in  the  simulation 
software.  In  the  case  of  the  example  used  in  this  report,  the  first-in,  first- 
out  discipline  was  used,  and  the  corresponding  time  that  a CT  waited  for  the 
channel  was  averaged  over  all  CT's  generated  during  the  simulation  run  as  well 
as  over  only  those  CT's  that  experienced  some  delay.  These  average  figures  are 
recorded  in  table  1 together  with  the  total  number  of  CT's  generated  and  the 
number  of  those  that  did  not  encounter  any  delay.  In  addition  to  waiting  times, 
second-by-second  observations  of  the  queuing  state  can  be  used  to  derive  the 
number  of  CT's  waiting  for  the  channel  averaged  over  the  t th  minute  of  si  ula- 
ted  time  and  the  entire  sample  period  of  2 hours.  The  120-minute  average  are 
plotted  in  figure  8,  and  the  2-hour  average  is  provided  in  table  1 along  with 
the  maximum  number  of  CT's  in  queue  in  any  second.  Thus,  the  average  level  of 
the  graph  in  figure  8 is  .426,  and  the  maximum  of  the  queue  lengths  observed 
during  successive  seconds  of  simulated  time  is  6. 

APPLICATIONS  OF  THE  MODEL . 

The  model  has  been  used  to  evaluate  the  effect  of  a proposed  FCC  regulation  on 
air  traffic  control  air/ground  communications  (references  4 and  5).  The 
regulation  required  transmitter  identification  in  the  form  of  coded  tone  bursts 
at  the  beginning  and  end  of  pilot  initiated  TR's  from  privately  owned  aircraft. 

The  duration  of  each  tone  was  anticipated  to  be  between  a few  tenths  of  a second 
and  1.5  seconds.  This,  of  course,  is  equivalent  to  increasing  the  length  of 
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TABLE  1.  GPSS  SIMULATION  MODEL  FOR  ATC  VERBAL  COMMUNICATIONS  SYSTEM 


FIGURE  7.  SIXTY-SECOND  AVERAGE  CHANNEL  UTILIZATION 


LOADINGS 


pilot-initiated  TR's.  These  increases  were  incorpororated  into  the  model.  In 
particular,  simulation  experiments  were  conducted  for  tone  burst  lengths  of 
0,  0.5,  1.0,  and  1.5  seconds  and  the  effect  upon  the  basic  response  variables 
was  observed . 

In  another  application,  the  model  was  used  to  determine  the  effect  of  aircraft 
arrival  rate  on  sector  communications  (reference  4).  Obviously,  as  the  rate 
Increases,  the  Intensity  of  communications  traffic  and  the  number  and  length 
of  communication  delays  will  eventually  become  intolerable.  This  is  reflected 
by  the  model  in  terms  of  the  channel  utilization  and  the  aircraft  queuing  state. 
The  problem  is  to  ascertain  that  rate  or  range  of  rates  beyond  which  adequate 
air  traffic  control  service  cannot  be  maintained  because  of  demands  placed 
upon  the  communication  system.  In  this  sense,  sector  capacity  in  terms  of 
aircraft  .rrlval  rates  can  be  established  and  compared  with  rates  expected  in 
future  years  on  the  basis  of  air  traffic  forecasts. 

Other  applications  of  the  model  can  be  found  in  reference  4.  It  suffices  to 
say  here  that  the  ATC  communication  model  is  a highly  versatile  analytical 
tool.  In  its  present  form,  the  model  can  be  used  to  evaluate  the  effect  of 
changes  in  any  of  the  input  variables  on  sector  communications.  Moreover,  by 
using  historical  records  obtained  from  operating  sectors,  the  model  may  be 
adapted  to  situations  not  covered  by  its  present  software  realization.  For 
example,  the  model  does  not  presently  distinguish  between  pilot  and  ground 
initiated  TR's.  However,  as  in  the  case  of  the  tone  burst  study,  it  is  possible 
to  determine  from  historical  records  (e.g.  see  reference  1)  the  proportion  of 
TR's  initiated  by  the  pilot.  Thereafter,  it  is  a relatively  simple  matter  to 
modify  existing  software  so  that  the  model  does  distinguish  between  pilot  and 
ground-initiated  TR's. 


CONCLUSIONS 


An  ATC  communications  simulation  model  has  been  realized  in  the  form  of  com- 
puter software.  It  is  a versatile  analytic  tool  for  studying  second-by-second 
behavior  of  air/ground  voice  communications  in  a specific  control  sector,  say 
New  York  LC  sector  510,  or  a class  of  control  sectors,  say  a typical  HI  sector 
in  the  Houston  Air  Route  Traffic  Control  Center.  Input  parameters  include  air- 
craft arrival  rate  and  distributions  of  TR  lengths,  number  of  CT's  per  aircraft, 
and  number  of  TR's  per  CT.  The  principal  outputs  include  second-by-second 
observations  of  the  number  of  aircraft  in  sector,  channel  utilization  (busy  or 
idle) , and  aircraft  queuing  state  (waiting  for  an  already  occupied  channel  or 
otherwise.)  The  model  is  presently  available  as  usable  software  in  computer 
facilities  at  NAFEC. 
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APPENDIX 
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PFRCFNT  OF  /FRO  ENTRIES  = f01,6/XXX.X# 

AVFPAGF  NUMRFR  OF  AIRCRAFT  IA  OllFIIF  = 801 , 3/XXX.XXX8 

MAXIMUM  NUMRFR  OF  AIRCRAFT  IN  OllFIIF  = 8 ni,?/XXX8 


•*xxxxxxxxxxx»xacexecs«cxxxxx  FND  xxxxxx*xxxxxxxxxxxxx«xx*cxxxxxxxxxxxxxxxxx 


FNO 
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^ I IB  B PUT  INF  FXPnM  F X , I X , H ,1?  , I ?, 14  I 
R = R4Mnil(  !1  ) 

X = ( -FX ) *A LOCI R ) 

IXrX+O.S 

RFTIIBN 

F«n 


FMF'C  T IHM  PAMTIH  ] X ) 

JYr IXSAFS30 
1 F ( J Y ) 5 » 4 , 4 

3 I y=  ]v+?i  474P3A47+1 

iS  RAMP\|=1Y 

RA  MO|  1=  R 4 FilTUS  .4  AS  A4]  3F-R 

|X=1V 

RFTIION 

FFtn 


SIW  PH  IT  IMF  ftIBNPl  (X  ,P,N'X,  11  , I ? » I ? ) 

IMPLICIT  RFALSPID) 

R FA  L K 
OP  = P 
DK  = K 

R = RAMOII(  II  ) 

PP  = R 
I = n 
x = n .0 

nr.  iiM=np««nx 
n y = nr.  iim 

IF  ( ncii“  .it.  hp  ) r,n  tp  io 

wx  = ( x+o.f  ) 

R F Tt |R  N 
10  1=1+) 

X=  1 .0#I 

OX=X 

nY=PY*((nK  + nx-i  .onm/nxi 
nr  i im=  nc  i im+  o y * ( 1 . o on-  np  ) » *n  x 
if  i rc  »IM  .IT.  OR)  r.p  to  io 
NX=( X+O.S) 

RFTMRN 
FMO 
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I 

St'BPPIIT  IMF  f APMIM  ,A?  ,N  ,Y  , M , 12  ) 
C"=*.0 

j CLO.l 

XN=A)*A?*M 

if  ( a?  ,i f.  o . i r,n  m i 

j SO=(CII-X“  1/2.F7SR 

r,n  Tn  2 

1 S0=l  XM-r.L  1/2.S7S8 

2 X*RMPRM(  I 1 ) 

XX  = ( X*sr>)+ X* 

I f ( xx  .if  . o. ) xx=n. 

y=Fyp ( XX ) 

R FTIIPM 
F^n 


FUNCTION  RMTRMI  IX  ) 
S"M  = 0 .0 
nn  f 1=1,1? 

; SHMsStlM+RAMPtll  IX  ) 

RMriRF=SI'M-f  .0 

j RFTIIRN 

i fMP 


$ I IR  801  IT  IMF  MSRr.t  PI  , P?  ,M  ,MT  I MF  ,1  1 , 12) 
CALI  SIJPf.AM(  Pi  ,P2  ,X  , I 1 , 12  ) 

NT  JNFaX  + O.S 

IF  ( NT  IMF  .IT . I ) NT IMF  = l 

RFTIIRN 

fMD 


I 

< 


o n o o o o r>  o o o nnnnonoon 


Surpoiit  IMF  Simr.AMf  x ,a  ,x  , I ?,  14  ) 


RFAl  K 

C 

C TH]S  SI  iRPpi  IT  I HP  GFNFRATFS  RANDOM  VAR  I A T FS  FR  X A G A MM  A DSN.  WITH 

C GFNFPAL  FORM 

C X K-l 

C A X F XPI  -AX ) 

FIX)  = 

GAMMA ( K ) 

WHFN  X IS  AM  IMTFT.FR,  T HF  GAMMA  nSM.  IS  FDIIIVAI.FNT  TO  THF  FRLAMG 

n f m . , whit.h  abjffs  af  a sir  of  k exponential  vapiatfs  vith  fxpfttfd  valmf 

1/A.  THFRFFOPF,  THF  FRL  AMO  V A R I A T F X IS  EQUAL  TO  1/A  T1MFS  THF  LOG  OF  THF 

ppnniCT  of  x PAMnoR  vapiatfs  fro«  a iinifirm  in,])  dsn.. 

WHFN  X is  MPT  AM  [MTFT.FR,  AM  APPROXIMATF  TFT.HMjnijF  MOST  RF  IISFD. 
lft  K = M+P  WHFPF  M IS  THF  SMALLEST  imtft.fr  com  t a i mf  n in  K AMn  0 IS  THF 
REMAINDER.  SIA'TF  THF  FXPFTTFn  VALOF,  V AR  I A Mf.F  , AMn  T H I P n TFMTRAL  MOMENT 
OF  GAMMMA  VARIATFS  APF  LINEAR  FIIMT.TinNS  rF  X,  AN  APRRPXIMATF  TFT.  HM  1 01  IF 
FDR  T.FNFRATINT.  GAMMA  VARJATFS  WITH  PARAMETER  X ]S  TP  GFMFRATF  A MIXTURE 
OF  GAMMA  VAPIATFS,  C HPOS  I MG  m WITH  PROP.  (1-0)  AMD  M+l  WITH  PROP.  0. 

THF  A PPRPX  IMAT  )PM  jMppnv/FS  WITH  INCREASING  K. 

REF.  - ' C pm  PI  IT  F P SIMULATION  TECHNIQUES' 

MAVLOP,  RALIMTFY,  ROPQICK,  f.  CHI ) 

JOHN  WHEY  f.  SONS,  1 PGA  PP  . 87-90 

Ml  *X 
m?=mi ♦! 

0 = K- FL  OA  T (Ml  ) 

KX  *M] 

IF  (P  . FO . o.n)  GO  TO 
T * R A MOO I IT  I 
IF  (T  ,LF  . O)  KX  = M2 
10  TP*1.0 

oo  ?n  1=1 ,xx 

R = R A MOlll  14  ) 

20  tr=tr*r 

X*(-ALOG(TP  ) ) / A 
RETURN 
FMO 
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FUR  BOUT  IMP  T 1MFSI  1 VAt  (IF,  JCAVFF,  ISAVFH,  l FAC,  I STn.FFTP,  JOllF, 

* FOIIF  , 1100  , I TAR  ,F1  AR  , JIISF,  I IKCF  .FIISF  » IPAX  , I *AXO  , I MAXH,  JMAXPH, 

* FSA  VFL  , I W/>  XL  ,FNA  XRL  1 , 

1 MTFr.FR  s ? ] s AVFH  , II  nr  , I IlSF  , I MAX  RH 

rfal*R  f n* ip  ,FHSF 
R F A t * A FSTP.FRAVFl  , F F»A  XR  |_ 

OJMFA'S  I pm  ; VALIIF(  A ),  1 SAV/FFI2  ) , ISAVFHl  2 ) , IF  AT.  ( ? ) , [ S TOf  2 ) , F STP(  2 ) , 

J IFMIF(  2)  ,FOHF  ( ? ),  ILnr.l  ? ) , ITARI  ? ) ,FTAP(  2)  ,1  USF(  ? ) , lUSFFI  2 I ,FIISE(  2 ) , 

* 1 PA  x ( 2 ) , I MA  XR  I 2 ) , 1 *'A  XMI  ? ) , lMAXPHl 2 I ,F  SAV  FLI  2 1 , I M A X L I 2 ) , FM A X P L ( 2 1 
1MTFP-FR  Z«  POO)  ,KFYC  ? )/I  ,10, 12/ 

MI1Kin=  I FA  VFH(  1 ) 
ir  sf.PUlMM  MMMPFR 
I D s Q pw  Ml  IMP  FR 
I C = 1 
IfM-  1 

K = Fi 

*■ — 

L ■=  1 

n n 2P0  i jk  = i ,3 

M s K F V ( I JK  ) 

JK=K*(N-1  ) ♦ (_ 

OP  100  1 R =1  , POO 

J-l  J*AXH( JK 1 +2*1  IO'*( 1R-1 )♦(  IC-1 1 ) 1/2 

100  21 IR)=IMAXRM( j) 

wo  I TF  ( Ml  (MI  T , 10  1)  Z 

101  FPRMaTIAO]?) 

’00  COMTINUF 

‘ PFTIIPN 
FNO 


\ 

I 


I 


i 


I 


51  IB  POUT  I MF  PA5  52(J  VAl  IIF,  15AVFF  , I SAVFH,  I FAC  , I 5Tn,F5TO,  IOIIF, 

* FPIIF  , 1 1 nr,  , MAR  , F1AR  , IIISF,  I U5FF  ,F(ISF  t I MAX  , I MAX  R , I MA  X H , I M * XR  H , 

* F5AVFL  , IMAXL  ,FRAXRl  ) . I 

iNT ft.fr*?  t«, avfh , iLnr-, mse , ihaxbh 

R FA  t «R  F ni  IF  ,FI'5F 

RFAL*A  FSTO.FFAVFL .FRAXBl 

0 J«FM5 ION  ? VAU'FI  A ) , I4A  VFFI2  ) , ISAVFHt  2 ) , 1FACI  2)  , ISm<  2 ),FSTn<  2)  , 
.*  iniiFt  ?)  ,FOHF  ( 2) , lllir.l  ■>  ) , | TAB  ( 2 ) ,FTAR(  2 ) , | US  F ( 2 I , I IlSFFt  2 I . FI  IS  6 ( 2 ) < 

* ] MA  X ( 2 ) , I MA  >P ( 2 ),1«AXH(?1  ,]MAXFH(?),FSAVFL(  ?)  ,|M*xi(2l,FMAXBL(  2) 
JMTFr.FR  ;i  1 200  ) ,M|  <J  ) ,KFYM|  9 ) / 3, 3,4,5  ,4,4,  7 ,R  ,9/ 

t MT  F OF  R KF  VM(  41/2,3,4, 5, 4,7, B, 9,11/ 

NUN ! T = I S A v F h I 2 ) 
ir  = i 
ir  m«i 
Kx  A 

L=  1 ^ 

nn  loo  1=1,9 

KMxKFYMI i t 

100  M(  n x i s a v F M ( KM » 

WP  ITF  I Mt IM I T , IP  1)  N 

101  F PP  RA  T ( 9 I 5 ) 

OP  200  I = 1 ,9 
KMxK  FYM(  I ) 

NP=M(  t) 
jkxk*(km-i t+L 
on  300  I R = J , M R 

J = ( JMAXWI  JK  ) + ?*(  ICM*(.IB-l.V*(  IC-llll/2 
300  2 ( I R I x imaxbh(.i  ) 

WOITFI  NIIMIT  ,102  ) ( 7 ( M ) ,M  = 1 ,MR  ) 

102  forma  T ( 1 A]  5 I 
700  CPNTINUF 

R F Tl  ION 

EMO  . _ 
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I 

. -A 


L 


SI  « urn  IT  INF  FXPOM  F X , IX  , II  ,1?  , I 14  ) 
r*rawo(i( n ) 
x<i - f x > *Ainr.<  r ) 

IX*X*0.S 

RFTIIBN 

F»m 


FUNCTION  ► T>l  MIX) 

I v*  IX  #.',SS3 p 
IF!  IY)5,A,S 
IY=lY4?l<.7Af>3AA7  + l 

HANOI'*  I Y . • jJ'  vjjiTt. 

R A NOl '*  R A NOI  i*  . 4 AS  AM  3F-9  ■ \ (j  I *** 

IX*  IV  » J*~m* 

RFTHON 

ewn 


SIWRMIT  IMF  StlRMPl  (K  ,P  ,MX,  II  , I ? 1 1 ? I 
IMPLICIT  R F A L * P ( 0 ) 

R FA  L K 

OP*P 

ok*k 

RxRAMDU!  II  I 
DR  * R 
l*n 
x*o.o 

DC  |IM=nP**DK 

DY-nciiM 

IF  i nciiw  .it.  hr  ) r,n  td  io 
MX* ( X+0.5  ) 

RFTiibn 
I- !♦  1 
X*  1 .0*1 

ox*x 

nv»nY*( ( ok+ox-i .onn )/nx ) 
oriiM*ncMM4  0Y*(  l .non- op  )**nx 
IF  IPCIIM  .LT.  OR)  r.n  TCI  10 
NX*( X+0.5) 

RFTIIBN 

fwo 


A-ll 


simpni it imf  papm(ai  ,a?,n,v  tn  1 12) 

c • <s . o 

XM=A]+A2*N 

ifia?  .IP.  o.)  r.n  to  i 

snrir.u-XM 

r.n  t n 2 

1 SP«( >M-CL 1/2.5758 

2 X*RNPRW(  I 1 ) 

x x * ( x*SP)+  xM 

IF(XX  .LF  . 0.  ) xx=n. 

Y=FXO(  XX  ) 

RFTIIRN 

pwn 


FUNCTION  RMPRMIIX) 
SlINsP  .0 

np  5 I«1 , i? 

5 SlIMsSl'M+RANnill  IX  ) 

. RMnRN=SHM-A  .0 
RFTIIRN 
FNO 


SI  I*  ROOT  IMF  MSSr.l  PUP?  ,N,MT  IMF,  I 1,12) 
CALL  Simr.AMiPi  ,P2,X,I1  ,12  ) 

nt  jMF»x  + n.s 

IF  (NT  INF  ,LT  . I I NT IMF  = 1 

RFTIIRN 

END 
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I 
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SI*  «n(IT  IMF  S(|RPAM<  X ,A  ,X  , I ? , 14  ) 


RFAl  K 

THJS  SI  IRP  PI  IT  I MF  PFNFRATFS  RAMPpM  VAR  I A T FS  FR  rM  A PAMMA  l)SN.  WITH 
PFNFPAL  FORM 

K K-l 

A X EXPI-AX) 

FIX)  = 


PAMMA  I K ) 

WHFN  X IS  am  IMTFT.fr,  TMF  GAMMA  ITSM.  IS  Fnnl  VAI.  FVJT  TP  THF  FRLAMp 
PSM , , WHICH  A P I S F S AS  A SI'H  PF  X F XPPNFN  1 I A L VAPIATFS  ||TH  FXPFCTFP  VALPF 
1/A.  THFPFFPRF,  THF  FPL  AMP  VAR  I A T F X IS  FIJI  I A L TP  l/A  T 1 M F S THF  LOG  PF  THF 
PPPPICT  PF  X PAMOPM  VAPIATFS  FPP«  A IIMJFIRM  IP,  1)  PSM.. 

WHFM  K IS  MPT  AM  I MTFPFR , AN  APPRflXI  MATF  T FC  HM|  PI  IF  MUST  R F IISFO. 

L FT  X - M+  P WHFBF  M JS  THF  SMALL  FS  T INTFPFR  CPM  T A 1 MR  P ]M  K AMP  0 IS  THF 
QFMA  I MflFP  . SIA'CF  THF  FXPFCTFP  V A L 1 1 F , V AR  I A MCF  , AMO  THIPP  C FNTRA  L M QM  F N T 
PF  PAMMMA  VAPIATFS  A P F L1AFAR  FlIMf.  T I PMS  fF  X,  AM  APPRPXIMATF  TFC.HNICIIIF 
FnR  PFNFRATJMP  PAMMA  VAPIATFS  WITH  PARAMFTFP  X IS  TP  C.FMFPATE  A MjxTlIHF 
PF  PAMMA  VAPIATFS,  CHOOSING  M WITH  PRP«.  (1-0)  A NO  M+l  WITH  PROP.  0. 

THF  APPROXIMATION  IMPPPVFS  WITH  IMCPFASING  K. 

REF.  - * C Pmpi  /T  F P S I M | il  A T I ON  TECHNIQUES* 

NAYLPP,  PAL  IMTFV,  RiiPPICK,  f.  C.HII 
JOHN  W I L F Y f.  SONS  , 1944  PP  . H7-90 


M1=K 

M?  = Ml  ♦ 1 

P = K-  FL  OA  T I M 1 ) 

KX  =M) 

if  ip  ,fo.  r.m  p n to  10 
t = r a mpiii  n i 

IF  ( T .IF.  O)  XX  = M? 

10  T R = 1 . 0 

r>n  ? o i = i ,xx 

R.RAMplll  14  ) 

20  TP=TP*R 

X* I -A  LOG! T P ) ) / A 

RETURN 

EMP 


nw 

> « « i .« w^riULL  LUr.J 
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Si  IP  POUT  IMP  T IMPS!  I VALUF  , I S A VF  F , I S AVF  H , I FAC  , I STO  ,F  5 TO  , JCU)  F , 

* FOIIF  , I LOO.  ,1  TAR.FTAR  , I MS  F , 1 1 IS  F=F  ,F(jSF  , I MAX  , I ma  XR  , I MA  X H , I MA  XR  H , 

* FSA VFL , 1 MA XL  .FRAXRl  1 ' ( 

IMTFr.FRe?  isavfh , uno,  IIISF , IMAXBH 

RFAL°R  FOIIF  , F MS  F 
R F A L * A FSTP.FSAVFL .FMAXRL 

0 IMF  MS  I PM  I VALMFI 6 1 , I SAVFFI 2 ) , ISAVFHI 2 1 , 1FACI  2 ) , ISTOI 2) ,FSTPI 2) , 

» IOHF(  2)  ,FOUF(2  1,  ILPOI 2 ) , ITARI2 1 , FTAfM  2) ,1 USFI 2 ) , IUSFFI 2) ,FIISE( 21 , 
•IMAXI2)  , IMA  X«( 2 ) , IMAXMI 2 ) , IMAXRHI 2 I tFSAVFLI  2 » ,IMAXL( 2) , FMAXPU  2) 

IMTFr-fR  Z(  AOO  ) ,K  FYI  ? )/l  , 10,  12/ 

Ml  IM  I T = I S A VF  HI  1 ) 
ir.=cninMM  mumrfr 
] R = R PW  MIIMRFR 
IC  = 1 
1TM=1 
K = 6 

1 = 1 

r»n  200  I jk  = 1 ,3 
M = K F V(  I JK  ) 

JK=K*(M-] )+L 
OP  100  I R = 1 ,600 

J=(  IMAXHI  JK  ) *2e(  ICM*(  JR-l  )♦(  IC-1  m/2 

100  z ( IR ) = IMA XRHI J ) 

WR  1 T F ( NIIMI  T , 10  1 1 2 

101  FOR  MA  T I 40 12) 

?00  CPMTINIJE 

RFTIIPN 


